
function[M2N]=G_matrix_2XN(xi,xb,yb,N,mu,epsilon)
M2N=zeros([2,2*N]);    
for j=1:N
        xj=[xb(j), yb(j)];
        Rij=(xi-xj);
        Rij2=sum(Rij.^2);
        rhoij2=Rij2+epsilon^2;
        rhoij=sqrt(rhoij2);
        Aij=1/(4*pi*mu)*(log(rhoij+epsilon)-epsilon*(rhoij+2*epsilon)...
                       /(rhoij*(rhoij+epsilon)));
        Bij=(rhoij+2*epsilon)/(4*pi*mu*rhoij*(rhoij+epsilon)^2);
        Gxx=-Aij+Bij*Rij(1)^2;
        Gxy=Bij*Rij(1)*Rij(2);
        Gyx=Gxy;
        Gyy=-Aij+Bij*Rij(2)^2;
        M2N(1,(j-1)*2+1)=Gxx;
        M2N(2,(j-1)*2+1)=Gxy;
        M2N(1,(j-1)*2+2)=Gyx;
        M2N(2,(j-1)*2+2)=Gyy;            
    end











